LIBRARY 

SBCTKW 

JE  sCHOCft 
MONTEREY.  CAUFOBNIA    SOB* 


NPS-69-78-012 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


A  COMPUTER  SUBROUTINE  FOR  STRESS  ANALYSIS 

OF 

ROTATING,  HEATED  DISKS 

by 

John  E.  Brock 

Robert  E.  Brown 

May  1978 

Approved  for  public  release;  distribution  unlimited. 
Prepared  for: 


FEDDOCS 

D208.i4/2:NPS-69-78-0i2    Chief  of  Naval    Research 

Arlington,   Virginia     22217 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  Tyler  Dedman  J.  R.  Bor sting 

Superintendent  Provost 


A  COMPUTER  SUBROUTINE  FOR  STRESS  ANALYSIS 
OF  ROTATING,  HEATED  DISKS 


This  report  gives  listing  and  instructions  for  using  a 
digital  computer  subroutine  for  finding  stress  distri- 
bution in  a  thin  rotating  disk  with  nonuniform  heating; 
the  problem  is  axisymmetric.  An  iterative  method  is  used. 
Theoretical  background  is  given. 


UNCTASSTFTF.D 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1.     REPORT  NUMBER 

NPS- 69-78-012 

2.  GOVT  ACCESSION  NO. 

3.     RECIPIENT'S  CATALOG  NUMBER 

4.     TITLE  (and  Subtitle) 

A  COMPUTER  SUBROUTINE  FOR  STRESS  ANALYSIS  OF 
ROTATING,   HEATED  DISKS 

5.     TYPE  OF  REPORT  &  PERIOD  COVERED 

6.     PERFORMING  ORG.   REPORT  NUMBER 

7.     AUTHORfs,) 

John  E.  Brock 
Robert  E.   Brown 

8.     CONTRACT  OR  GRANT  NUMBERfs) 

9.     PERFORMING  ORGANIZATION    NAME   AND   ADDRESS 

Professor  John  E.  Brock  (Code  69Bc) 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 

10.     PROGRAM   ELEMENT,  PROJECT,   TASK 
AREA  &  WORK  UNIT  NUMBERS 

611 52N,   RR000-01-01 
N0001478WR80023 

11.     CONTROLLING  OFFICE  NAME   AND  ADDRESS 

Naval    Postgraduate  School 
Monterey,   California     93940 

12.     REPORT  DATE 

May  1978 

13.     NUMBER  OF  PAGES 

14.     MONITORING   AGENCY  NAME  &    ADORESSf//  different  from  Controlling  Office) 

Chief  of  Naval   Research 
Arlington,   Virginia     22217 

15.     SECURITY  CLASS,  (of  this  report) 

15a.     DECLASSIFI  CATION/ DOWN  GRADING 
SCHEDULE 

16.     DISTRIBUTION   ST  ATEMEN  T  (of  this  Report) 

Approved  for  public  release;  distribution  unlimited. 

17.     DISTRIBUTION  STATEMENT  (of  the  abstract  entered  In  Block  20,  It  different  from  Report) 

18.     SUPPLEMENTARY  NOTES 

19.     KEY  WORDS  (Continue  on  reverse  side  It  necessary  and  Identify  by  block  number) 

Stress  analysis 

Rotating  disks,  Heated  Disks 

Axisymmetric  Elasticity,  Axisymmetric  Disks,  Elastic  Disks 

20.     ABSTRACT  (Continue  on  reverse  side  It  necessary  and  Identity  by  block  number) 

This  report  gives  listing  and  instructions  for  using  a  digital  computer 
subroutine  for  finding  stress  distribution  in  a  thin  rotating  disk  with 
nonuniform  heating;   the  problem  is  axisymmetric.     An  iterative  method  is 
used.     Theoretical  background  is  given. 

DD  ,  ™»n  1473 


EDITION  OF   1  NOV  65  IS  OBSOLETE 
S/N   0102-014- 6601   | 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


TABLE  OF  CONTENTS 

Introduction 1 

Fundamental  Analysis 2 

Power  Law  of  Thickness  Variation 3 

Exponential  Law  of  Thickness  Variation 4 

General  Law  of  Thickness  Variation 5 

Computer  Implementation 8 

References . 9 

Appendix  A 10 

Appendix  B 13 

Appendix  C 16 

Initial  Distribution  List 20 


A  Computer  Subroutine  for  Stress 
Analysis  of  Rotating,  Heated  Disks 

by 
John  E.  Brock  and  Robert  E.  Brown 


Introduction 


Although,  as  is  indicated  by  the  title  hereof,  the  principal 
purpose  of  this  monograph  is  to  present  tested  and  proved  digital  com- 
puter software  for  the  analysis  of  stress  in  a  spinning  axisymmetric 
disk  having  a  radially  variable  thermal  strain  field,  the  opportunity 
is  also  taken  of  developing  the  theory  and  presenting  some  analytic 
solutions. 

The  method  developed  herein  for  computer  analysis  of  disks 
having  a  general  law  of  thickness  variation  was  suggested  by  the  al- 
gorithm contained  in  reference  2  and  it  appears  to  have  advantages 
over  such  procedures  as  that  of  M.  Donath,  reference  3,  x^hich  has  been 
widely  circulated  in  a  book  by  S.  Timoshenko,  reference  5. 

In  what  follows  we  immediately  obtain  a  second  order  linear 
differential  equation  with  dependent  variable  u,  the  radial  deformation, 
and  r,  the  radius.  Analytic  treatment  is  given  for  two  particular  laws 
of  thickness  variation.  For  the  general  case  of  thickness  variation, 
the  equation  is  recast  as  a  second  order  linear  differential  equation 
in  which  the  dependent  variable  is  radial  stress,  a   .  However,  for 
numerical  treatment  an  alternate  form  is  more  useful  and  direct  and 
this  forms  the  basis  of  the  digital  computer  software  which  is  given 


and  illustrated  in  the  appendices  hereto. 

Fundamental  Analysis 

We  presume  that  the  disk  is  thin  enough  and  that  the  thick- 
ness varies  slowly  enough  with  respect  to  radius  that  we  are  just- 
ified in  neglecting  all  stress  components  excepting  only  the  radial 
stress  a  and  the  circumferential  stress  ofi .  Material  properties 
E,  Young's  modulus,  and  v,  Poisson's  ratio,  are  presumed  to  be  in- 
deed constant.  The  thermal  strain  field,  aT,  and  the  density  y 
may  be  specified  functions  of  radius. 

Two  types  of  problem  are  considered: 

1.  Annular  disk,  0  <  a  ^  r  =  b,  with  a   (a)  and  a   (b) 

r        r 

being  specified. 

2.  Solid  disk,  0  i  r  =  b,  with  a  (b)  being  specified. 
We  also  use  the  symbols  t  =  t(r)  for  disk  thickness  and  ca  for 
angular  velocity.  Other  simplifying  notation  will  be  introduced 
later  on. 

Consideration  of  radial  equilibrium  leads  without  difficulty 

to  the  equation 

1  d(rta  )         22  f-i  \ 

r- n_  -  Of,   -  yarrz  (1) 

t      dr      e 
The  thermoelastic  constitutive  equations  are 

EeQ  =  aa  -   va  +  EaT;  Ee  =  a  -  vaQ  +  EaT        (2a, b 
9    9     r     '    r    r    9 

where  the  strain  components  are 

eQ  =  u/r;   er  =  du/dr  (3a,b 

Strain  compatability  leads  to  the  equation 

r  ^(Eu/r)  =  -(l+v)(a9-ar)  (4) 

These  equations  may  easily  be  combined  into  the  differential  equation 
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u"  +  u'/r  -  u/r2  +  (t'/t)[u'  +  vu/r  -  (l+v)aT] 

=  (l+v)(aT)'  -  (l-v2)yw2r/E        (5) 
where  primes  denote  differentiation  with  respect  to  r. 

Two  particular  laws  of  thickness  variation  permit  siirrole 
analytic  treatment. 

Power  Law  of  Thickness  Variation 


If 


t  =  t0(r/r0)n  (6) 


so  that 


(f/t)  =  n/r  (7) 

then  equation  5  becomes 

u"  +  (l+vn)uVr  +  (vn-l)u/r2  =  8'  +  nB/r  -  kr         (8) 

where 

3  =  (l+v)ctT,     k  =  (1-v2)yoj2/E  (9,10) 

Equation  8  may  be  rewritten  as 

d  r  1-m  d  ,    (m+n)/2  Nn    (2+n-m)/2  ,_,    n/       ,    N   ,,,v 
^[r    ^r  u)]  =  r'        (8'  +  nS/r  -kr)  (11) 

where 

m  =  ±/(n2-4  n+4)  =  ±/[(n-2)2+Ml-  )n]  (12) 

and  either  the  positive  or  the  negative  sign  may  be  used  Equation  11 
may  be  proved  simply  by  performing  the  indicated  operations  and  com- 
paring with  equation  8. 

The  quantity  on  the  right  in  equation  11  is  well  defined  so 
that  the  solution  of  the  differential  equation  may  be  obtained  simDly 
by  integration,  multiplication  by  r  ~  ,  and  another  integration.  Two 
constants  of  integration  are  introduced.  For  the  solid  disk  (case  2) 
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u(0)  =  0  gives  one  condition  and  the  second  comes  from  satisfying 
the  given  value  of  a  (b).  For  the  annular  disk  (case  1)  satisfying 
the  given  conditions  a  (a)  and  a   (b)  permits  evaluating  the  con- 
stants. 

An  example  with  n  =  -  0.42  is  given  in  Appendix  C.  Note 
that  the  case  n  =  0  corresponds  to  a  disk  of  uniform  thickness. 

Exponential  Law  of  Thickness  Variation 

If 

t  =  t0  exp(-mr2)  (13) 

where  m  is  a  constant  of  appropriate  dimensionality,  then 

(t'/t)  =  -2rm  (14) 

and  equation  5  becomes 

u"  +  (1/r  -  2rm)u'  -  (1/r2  +  2vm)  =  -2rm8  +  B'  -  kr     (15) 

If  additionally  we  assume  that  $ '  =  0  (which  makes  the  thermal 

strain  field  constant  —  a  triviality)  the  solution  is  simply 

u  =  r(fHk/2m)/(l+v)  (16) 

In  this  case  we  can  easily  find 

o„  =  aa  =   yw2/2m  (17) 

r    o 

which  is  independent  of  r.  Thus,  if  an  allowable  normal  stress  a. 
is  specified  and  if  blade  or  bucket  loading  at  r  =  b  is  w  (oounds, 
say)  per  unit  circumference,  then  a  disk  having  thickness 

t  =  (w/aA)  exp[(b2-  r2)(yw2/2aA)]  (18) 

will  be  such  that  a  =  aa  =   a. .  If  the  failure  criterion  is  the 

r    8    A 

maximum  shearing  stress  criterion  (Tresca's  condition),  it  is  clear 
that  this  disk  is  optimal  in  the  sense  of  having  least  volume  and 
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thus  having  least  weight.  If  the  failure  criterion  were  that  of 
von  Mises,  a  slightly  lighter  disk  would  suffice. 

General  Law  of  Thickness  Variation 

Equations  3a  and  2a  give 

Eu  =  rEaT  +  r(aQ-var)  (19) 

and  by  differentiation 

Eu'  =  EaT  +  rE(aT)'  +  (oa-va   )  +  r(a  '-va  ')      (20) 

y   r      y    r 

From  3b  and  2b  we  also  have 

Eu'  =  EaT  +  a  -  van  (21) 

r    9 

Subtracting  21  from  20  and  rearranging  gives 

E(aT)»  +  (l+v)(cQ-ar)/r  +  cjq  »  -  va^   =  0         (22) 
Equation  1  may  be  rewritten  as 

aa-a     =  ra  •  +  rva  +  Yw2r2  (23) 

8  r    r      r   ■ 

where  we  have  written 

v  =  t'/t  (24) 

for  convenience.  Differentiating  23  we  get 

a  '  =  2a  »  +  ra  "  +  va„  +  rva  '  +  rv»a^  +  2yw2r     (25) 
8     r     r     r     r       r 

and  substituting  23  and  25  into  22  gives 

ra  "  +  (3+rv)a  '  +  [(2+v)v+rv»]a  +  (3+v)yco2r  +  E(aT) '  =  0  (26) 
r         r  r 

This  is  a  single  differential  equation  with  dependent  var- 
iable cr  and  can  be  dealt  with  by  standard  numerical  methods.  The 
conditions  for  evaluating  the  constants  of  integration  have  been 
mentioned  earlier. 

However,  the  preceding  procedure  is  not  particularly  satis- 
factory. For  one  thing,  the  solid  disk  case  for  which  r  can  become 
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zero  encounters  numerical  difficulties  unless  special  treatment  is 
employed  to  evade  them.  More  importantly,  however,  if  T  is  given 
by  a  graph  or  a  numerical  table,  determination  of  (aT)'  may  involve 
numerical  difficulties  which,  ultimately,  are  due  to  our  having  per- 
formed the  differentiation  to  arrive  at  equation  20.  Accordingly, 
an  alternate  procedure  which  adheres  more  closely  to  the  fundamental 
mechanics  of  the  problem  is  described  in  what  follows. 

We  consider  the  annular  disk  first  and  we  represent  the  un- 
known stress  difference  on-<J     in  the  form 

8  r 

a9~°r  =  (A+Bn)r  (27) 

where  A  and  B  are  unknown  constants  and  ri  is  an  unknown  function 

normalized  so  that 

n(a)  =  0,   n(b)  =  1  (28a,b) 

Initially  we  make  an  assumption  for  n  taking  a  linear  var- 
iation in  the  absence  of  better  information.  By  the  use  of  some  of 
the  preceding  equations  we  will  be  able  to  construct  an  improved 
form  for  n  and  will  iterate  until  there  is  satisfactory  convergence. 

Let 

z  =  Eu/r,  w  =  EaT,   8  =  yw2  (29,30,31) 

noting  that  w  and  3  have  different  meanings  here  than  when  they 

were  used  earlier.  We  can  recast  equation  1  as 

d(ta  )/dr  =  (aa-a   )t/r  -  yco2rt  =  t(A+Bn-ycu2r)        (32) 
r       a  r 

Before  performing  the  indicated  integration  we  introduce  two 
convenient  notational  devices,  viz. 

p...  =  I  ...  dr;     *...  =  [...L=h  (33,34) 

Thus,  from  equation  32  we  have 
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to  =  (to)  +  Apt  +  Bpnt  -  pBrt  (35) 

r     r  a 

From  equation  4  and  equation  27  we  have 

dz/dr  =  -(l+v)(A+Bn)  (36) 


so  that 


Equation  2a  is 


so  that 


z  =  (z)  -  U+v)[A(r-a)  +  Bpn]  (37) 


z  =  w  +  (a«-o  )  +  (l-v)a  =  w  +  (A+Bn)r  +  (l-v)a   (38) 

o  V  V  V 


(z)  =  (w)a  +  Aa  +  (l-v)(c  )  (39) 

a     a  r  a 

(z)b  =  (w)b  +  (A+B)b  +  (l-v)(ar)b  (40) 

Evaluating  38  at  r  =  b  and  using  39  and  40  gives  one  equation 
involving  the  unknowns  A  and  B.  Evaluating  35  at  r  =  b  gives  a 
second.  These  equations  can  be  arranged  as 


*(pt)       *(pnt) 
(2+v)(b-a)    b+(l+v)*(pn) 


*(p8rt)  +  (tor)b  -  (tor)a   ] 

(W)  (W)  +(1_v)[(a  )  (a  )  ]j  W 


so  that  one  can  easily  solve  for  A  and  B.  Then  a     is  obtained  from 

35,  z  is  obtained  from  37  and  39,  and  (0^-0   )  and  ag  are  obtained 

from  38.  Then  a  new  function  n  is  calculated  from 

n  =  [ (aQ-ar )/r-A]/B  (42) 

Using  the  new  r\   the  process  is  iterated,  convergence  being 

monitored  by  examination  of  the  sequence  of  values  of  A  and  B  that 

are  calculated.  When  convergence  is  satisfactory,  the  desired 

functions  a  and  an   are  at  hand, 
r     8 

The  situation  is  simpler  for  the  solid  disk,  case  2.   (a  ) 
is  not  given  but  conditions  of  continuity  require  that  (o^-a   )/r 
vanish  at  r  =  0 .   Thus  A  is  zero  and  B  can  be  obtained  from 
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(w)0  -  (w)b  +  c[l-(t)0/(t)b](tor)b  +  c*(p0rt) 

B  =  b  +  (l+v)»(pn)  +  c*(pnt) (43) 

where 

c  =  (l-v)/(t)0  (44) 

Otherwise  the  procedure  is  as  for  case  1. 

Computer  Implementation 

The  theory  embodied  in  equations  27  through  44  and  the 
associated  procedure  has  been  programmed  for  digital  computer 
using  the  FORTRAN  language.  An  initial  programming  based  directly 
on  the  preceding  equations  and  written  in  January  1978  by  the 
junior  author  hereof  has  been  supplanted  by  a  newer  programming 
which  is  somewhat  more  compact  due  to  the  employment  therein  of 
ancillary  subroutines  developed  for  use  in  another  problem  (the 
lateral  buckling  of  elastic  beams)  on  which  we  are  working.  This 
program,  actually  a  subroutine  called  RODISK  is  listed  in  Appen- 
dix A  hereof.  This  listing  itself  contains  comments  which  ade- 
quately explain  the  construction  of  a  MAIN  program  which  supplies 
necessary  input  information  and  which  invokes  RODISK.  Appendix 
B  lists  the  ancillary  subroutines.  In  each  case  comments  indicate 
the  purpose  and  employment  of  the  subroutine.  These  may  prove 
useful  in  constructing  the  MAIN  program.  For  this  reason,  the 
ancillary  subroutines  DUFV  and  PRTV  are  given  even  through  they 
are  not  called  by  RODISK. 

In  the  computer  implementation,  the  various  functions   of 
r  which  appear  in  the  theory  are  represented  by  vectors   the  ele- 
ments of  which  are  function  values  at  equally  spaced  values  of  r. 
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The  ancillary  subroutines  manipulate  these  vectors.  All  of  these 
are  obvious  except  possibly  IMTV  which  performs  an  integration  by 
use  of  Milne's  formulas,  cf.  reference  k. 

In  the  subroutine  RODISK  there  is  a  slight  departure  from 
the  theory  as  given  herein.  As  a  first  step,  all  quantities  and 
functions  were  "dedimensionalized"  but  otherwise  the  method  is 
just  as  described  above.  Somewhat  finer  details  of  what  was 
actually  done  may  be  found  in  reference  1. 

Appendix  C  contains  some  examples  and  remarks  concerning 
them. 
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Appendix  A 
Listing  of  subroutine  RODISK 


(The  listing  on  this  page,  page  10,   is  of  the  comments  which 


provide  instructions  for  the  use  of  RODISK. 
on  the  following  two  pages . ) 


Commands  appear' 


SUBROUTINE    R 
THI S    IS    A    SU 
FEFENTIAL    ST 
ROTATING    AT 
ABOUT    ITS    AX 
DISTRIBUTION 
B'E    TREATED: 
TYPE    1. 
OUTSIDF 
INNER    R 
RACIUS 
TYPE    2. 
RACIAL 
THE    USER    MUS 
RODISK    AFTER 
(1)     N,     INTE 
INTO    WH 
DIVIDED 
D I  MENS  I 
BRAD 
ARAD       ( 
SRB 

SRA       (N 
TEEBEE, 
P0IS»    P 
KP(1)     = 
KP(2)  , 
OUTPUT. 
EVERY    F 
.  .  .  ,    9  6 
KP(3) , 
TO    BE    P 
ENGINEE 
KP (4). 
IF    KP(4 
PRINTED 
VECTORS 
VECTOR 
THICKNESS    PF 
RADII    STARTI 

VECTOR 
TIMES     (BRAD- 
GAMMA    DCES    N 
CONSTANT. 

VECTOR 
WHERE  (E)  IS 
LINEAR  THERM 
THE  MAI 
IMPLICIT  REA 
INTEGER  KP(4 
COMMON  X,N 
COMMON  /ONE/ 
FOLLOWI 
ANCILLARY  SU 
ON  THE  VECTO 
FROM  THE  LIS 
PROGRAM.  TW 
AVAILABLE  IN 
SUBROUTINE  P 
PRIV    WHICH    P 


ODISK.       JOHN    E.     BRCCK.,     1     MAY    1978 
BROUTINE    FOR    DETERMINING    RADIAL    AND    CIRCUM- 
RESSES    IN    AN    AXISYMMFTRIC    THIN    ELASTIC    DISK, 
ANGULAR    VELOCITY    CMEGA    (RADIANS    PER    SECOND) 
IS    OF    SYMMETRY    AND    HAVING    AN    AXISYMMETRIC 
OF    THERMAL    STRAIN.       TWO    TYPES    OF     PROBLEM    MAY 


(2) 
(3  ) 
(4) 
(5) 
(6  ) 
(7) 
(8  ) 
(9) 


(  10) 


(  11) 


(  12  ) 


ANNULAR  DISK  OF  INSIDE  RADIUS  ARAD  AND 
RADTUS  BRAD.   THE  RADIAL  STRESS  IS  SRA  AT  THE 

AHIUS  AND  SRB  AT  THE  OUTER  RADIUS.   THE  INSIDE 

MUST  PE  GREATER  THAN  ZERO. 

SOLID  DISK  OF  OUTSIDE  RADIUS  BRAD  AND  WITH 

STRESS  SRB  AT  THE  OUTSIDE  RADIUS. 

T  PROVIDE  A  MAIN  PROGRAM  WHICH  CALLS  SUBROUTINE 
IT  HAS  SUPPLIED  THE  FOLLOWING  INFORMATION. 

GER.   (iM-1)  IS  THE  NUMBER  OF  EQUAL  SUBDIVISION 

ICH  THE  ANNULAR  RACIUS  (BRAC  MINUS  ARAD)  IS 
FOR  COMPUTATIONAL  PURPOSES.   THE  PRESENT 

ONING  CAN  ACCOM MO CAT E  N  NOT  LARGER  THAN  10  1. 

NOT  NECESSARY  FOR  PROBLEMS  OF  TYPE  2.) 

OT  NECESSARY  FOR  PROBLEMS  OF  TYPE  2.) 

DISK  ^HICKNESS  AT  OUTSIDE  RADIUS 
OISSON'S  RATIO. 

1,2.   INTEGER  TO  DENOTE  PROBLEM  OF  TYPE  1,2. 
INTEGER  TO  PROVIDE  FOP  SKIPPING  WHILE  PRINTING 

FOR  EXAMPLE,  IF  N=101  AND  KP(2)=5,  ONLY 
IFTH  SET  OF  VALUES  WILL  BE  PRINTED:  1ST,  6TH, 
TH,  101ST. 

INTEGER  SPECIFYING  THE  NUMBER  OF  ITERATIONS 
ERFORMED.   USUALLY  KP(3)=10  IS  SUFFICIENT  FOR 
RING  ACCURACY. 

IF  KP(4)=0  ONLY  FINAL  ANSWERS  WILL  eE  PRINTED 
)=1  A  SEQUENCE  OF  ITERANT  VALUES  WILL  BE 
,  INDICATING  DEGREE  OF  CONVERGENCE. 

X(I  ,J) ,  1  =  1,2,3;  J=l,2,. ..  ,N. 
X(1,J)  CONTAINS  VALUES  OP  THE  RATIO  (LOCAL 

CISK)/(TEEBEE)  COMPUTED  AT  EQUALLY  SPACED 
NG  AT  THE  INSIDE  AND  ENDING  AT  THE  OUTSIDE. 
X(2,J)  CONTAINS  VALUES  OF  ( GAMMA) (OMEGA-SQUARE 
SQARED)  DIVIDED  BY  (SRB).   FOR  MOST  PROBLEMS 
OT  VARY  WITH  RADIUS  AND  THIS  QUANTITY  IS  A 


X(3,J)  CONTAINS  VALUES  OF 
YOUNG'S  MODULUS,  (ALPHA) 
AL  EXPANSION,  AND  (TEE)  IS 
N  PROGRAM  MUST  CONTAIN  THC 
L*8  (A-HtO-Z) 
) 


(F)  (ALPHA)  (TEE)/(SRB 

IS  THE  COEFFICIENT  0 
TEMPERATURE  CHANGE. 
STATEMENTS: 


ARAD, BRAD,SRA,SRB,TEEBEE,POIS,KP 
NG    SUBROUTINE    RODISK    THERE    ARE    SEVERAL 
BROUTINES    WHICH    PERFORM    VARIOUS    OPERATIONS 
RS    X(I,J).       THE    FUNCTION    OP    EACH     IS    OBVIOUS 
TING.       THEY    MAY    BE    USED    IN     THE    USER'S    MAIN 
0    OF    THESE    ANCILLARY    SUBROUTINES    WHICH    ARE 

THIS    PACKAGE    BLT    V»HICH    ARE    NOT    CALLEC    BY 
ODISK    ARE    DUPV    WHICH    DUPLICATES    A     VECTOR    AND 
RINTS    A    VECTOR. 
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s 
c 

c 
c 
c 
c 
c 
c 
c 

20   C 

c 
c 
c 

r 


B 

I 

7    F 

C 


USRO 

MPLI 

EAL* 

NTEG 

CNMC 

0MM0 

NE  =  1 

OIS  = 

ERO  = 

F(KP 

HC=A 

NM  =  N 

RITE 

ORMA 

0  5 
IM=I 
=  EIM 
(4,1 
(5,1 
(6,1 
T  =  R  = 
TA  =  X 
F(KP 
HE  P 
p.  AT  = 

1  =  (2 
ALL 
2=X( 
5=X( 
ALL 
ALL 
ALL 

6  =  X( 

ALL 

3=GN 

ALL 

ALL 

4  =  X( 

=  C1* 

=  (C5 

=  (C1 

F(KP 

QRMA 

ALL 

ALL 

ALL 

ALL 

=  CME 

ALL 

=  ETA 

ALL 

ALL 

J=A* 

ALL 

ALL 

ALL 

=  -(0 

ALL 

ALL 

=  POI 

ALL 

ALL 

ALL 

TER  = 

F(IT 

ALL 

ALL 

ALL 

0    TO 


)/ 

E0 


UTINE 
CIT    RE 
8    X(20 
Eft   KP( 

N    X,N 
N    /ONE 
.0  +  0 
2.D-1 
C.D+O 
1 1  ) .  EQ 
RAD/BR 

(6,2) 
T(//,l 

1=1, N 

-1 

/ENM 

)=RHO+ 

)=Y 

)=Y 

1 

(1,1 

(1). 

ROBLEM 

SRA/S^ 

.D+O+P 

INTV<  1 

7»N) 

3,1)-X 

MULVt  1 

i-ULVts 

INTV(9 
LC,N)+ 
INTV(6 
E+(ONE 
NULVt  1 
INTVd 
13,  N) 
C4-C2- 
*G4-C6 
*C6-C2 
(4). EG 
T(5X,I 
MULS(7 
VULSd 
ADDVt 1 
SUBVt 1 
-RHO 
MULSt  1 
*SRAT 
ADDSd 
DIVVI  1 
RHQ+X( 
MULSt  1 
I^ULS(5 
ADDVd 
NE-RHO 
MULSt  1 
AOOSd 
S-ONE 
MULSt? 
ADDV( 1 
SUBV(  1 
ITER+1 
ER.GT. 
OlVVt  1 
SUBStl 
CIVStI 
20 


RODI SK 

AL*3(A-H,0-Z) 
,101  ) 
4) 

E/ARAD,BRAD,S*A,SRB,TEE8EE,P0IS,KP 


.2)    ARAP=ZERO 
AD 

KPt  1  ) 

OX,' RODI  SK     PROBLEM    OF    TYPE     •  ,  I  1 ,  '  . ',  //) 


(ONE-PHO)*Y 


X  (  1 ,  M ) 

.2)     GO    TO    100 

IS    OF    TYPE    ONE:     ANMULAR    DISK 
p, 

0IS)*(ONE-RHO) 
,7) 

(3, N)-( ONE-POT S)*(OUE-S RAT  ) 

,2,8) 

,4,9  ) 

,10) 

(ONE-ETA*SRAT)/(QNE-RHO) 

,11) 

-RH0)*(0NE+P0ISJ*X<11,N) 

,6, 12) 

2,13) 


C3 

-03) 

*C5) 

.1) 

10,1 

tl4, 

3,15 

4,  15 

5,1  J 


/D 

/D 

!«RITE(6,7) 

P2E20.5) 

A) 

,S) 

,15) 

,16) 


ITER, A, 3 


6,16,S) 

6,16,S) 

6,  1,20) 

3, 1 )+(OME-PCIS )^SRAT 

1,11  ,3) 

,16, A) 

1,16,16) 

)*(0NE+POI S) 

6, 16, S) 

6,19,Z0) 

0,18,S) 
8,19,18) 
8,3,  19) 

KP(3)  )     GO    TO    200 
9,4, 17) 
7, 17, A) 
7,6, B) 
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100 


105 


200 


204 
235 


210 
211 


CI 
C2 
C5 
CA 
CA 
CA 
C6 
CA 
C3 
CA 
CA 
C4 
0  = 
SR 
8  = 
IF 
CA 
CA 
S  = 
CA 
CA 
S  = 
CA 
S  = 
CA 
S  = 
CA 
CA 

s  = 

CA 
CA 
IT 
IF 
CA 
X< 
CA 
GO 
CA 
CA 
CA 
CA 
S  = 
CA 
CA 
CA 
WP 
FO 
WP 
FO 
l'E 
KS 
DO 
J  = 
WR 
CO 
FO 
Ri 
EN 


=  P 

=X 
LL 
LL 
LL 
=  3 
LL 
=  0 
LL 
LL 
=  X 
CI 
AT 
(C 
(K 
LL 
LL 
ET 
LL 
LL 

-a 

LL 
PQ 
LL 
X( 
LL 
LL 
-B 
LL 
LL 
5R 
(I 
LL 
16 
LL 
T 
LL 
LL 
LL 
LL 
SR 
LL 
LL 
LL 
IT 
PM 
IT 
RM 
c# 

KI 

2 

1/ 

IT 
NT 
RM 
TU 

D 


IS-ONI 

A 

3tl)-X 

MULV(1 

MULV(7 

INTV(S 

E+X(9, 

INTVC6 

E+(ONE 

NULV(  1 

INTV(  1 

12,  N) 

C4-C2* 

(C5*C4 

-C6-C2 

(4).EQ 

V U  LS  ( 1 

SUBVd 

*SRAT 

AODS( 1 

0IVV(2 

(GNE+P 

VULS( 1 

S-ONE 

VULS(2 

,l)+(0 

ADDSd 

SUBVd 

(ONE+P 

NULSd 

ADDV(  1 

ITER+1 

ER.GT. 

MULS(4 

1)=0NE 

DIVVC1 

105 
NULS(4 
MULS(2 
MULSd 
ADDV( 8 
/BRAD* 
!*LILS<2 
NULS(3 
MULS( 1 
(6,204 
T(//// 
(6,205 
T(24X, 
LPHA.T 

=KP(2) 

0    1=1, 

SKIP 

(6,211 

NUE 

T(I10, 

N 


(3,N)+P0IS-CNF 

,2,7) 

,4,8) 

,9) 

N) 

,10) 

+POIS)*X( 10, N) 

,6,11) 

1,12) 

C3 

-C6*C3)/D 

*C5) /D 

.1)     kRITF(6,7> 

2,13,3) 

3,9, i^) 

4,20,S) 
0,1,20) 

01  S) 
0,15,S) 


0,18,S) 

NE-P0IS)*SRA7 

8,19,S) 

9,3,19) 

OIS) 

0,17,S) 

9,  17,19) 

KP(3)  )    GO    TG    200 
,16,6) 


ITER,SRAT,8 


9,16,6) 


,7 

0, 

9, 

,9 

*2 

,1 

,1 

,1 

) 

) 

) 

•R 

HE 

N, 

)J 

IP 


,BRAD) 
8,SRB) 
9,SRB) 
,9) 

0,S) 

1,SRB) 
5,TEEBEE) 


ADIUS*  ,1 IX,  •  THICKNESS  '  ,6X,  'GAMMA. OMEGA 
•  ,8X,' SIGMA.^AOI AL' ,7X,  •  S IGMA.C IRCUMF* 

KSKIP 

,X(7,I),X(15,I),X(1),I),X(11,I),X(£,I) 

6E20.5) 


.S< 

) 


,X( 
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Appendix  B 
Listing  of  ancillary  subroutines 


SUBROUTINE    ADDVCN1 »N2 fN3 ) 

C  ANCILLARY    SUBROUTINE 

C  ADC    TWO    VECTORS    TERM    3Y    TERM 

C  XIN3.I )=X(N1 ,I)+X(N2,I > 

R  E  AL  *  8    X  (  2  3  t  1  0 1  )  t  S 

COMMON    X.N 

DO    1    1  =  1, N 
1    X(N3,  I)=X(N1,I)+X(N2,  I) 

RETURN 

END 

SUBROUTINE     SUBV(N1,N2,N3) 
C  ANCILLARY    SUB^OUTIN^ 

C  SUBTRACT    TWO    VECTORS    T  ERM    BY    TERM 

C  X(N3fI)»X(NlfI)-X(N2fI) 

REAL*8    X(2J,131  )tS 

COMMON    X,N 

DO    L     1=1, N 
L    X(N3, I)=X(N1,I)-X(N2,I ) 

RETURN 

END 

SUBROUTINE    MULV (Nl , N 2, N3 ) 
C  ANCILLARY    SUBROUTINE 

C  MULTIPLY    TWO    VECTORS    TERM    BY    TERM 

C  X(N3,I)=X(N1,I )*X(N2.I ) 

REAL*8    X(20.101),S 

COMMON    XtN 

DC    1    1=1. N 
1    X(N3.I)=X(iMl,I)*X(N2.I  ) 

RETURN 

ENC 


C 


SUBROUTINE    DI VV  (Nl , N2 , N3 ) 
ANCILLARY    SUBPCUTINE 
f  DIVIDE    TWO    VECTORS    TERM    BY    TiRM 

C  X(N3.I)=X(N1,I)/X(N2,T) 

REAL*£    X(20,101).S 
COMMON    X.N 

DC    1     1=1. N  ,     „    T1 

1    X(N3,I)=X(N1,I)/X(N2,I) 

RETURN 
ENC 

SUBROUTINE    ADDS (Nl . N2 » S ) 
r  ANCILLARY    SUBROUTINE  ^,„r 

r  tor    &    SCALAR    TO    EACH    TERM    OF    A    VECTOR 

C  XIN2.I )=X(Ni ,1 )+5 

REAL*8    X(20,101).S 

COMMON    X.N 

DO    1    1=1, N 
1    X(N2.I)=X(N1.I)+S 

RETURN 

END 
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c 
c 
c 


SUBROUTINE    SUBS  (MU  N2  ,  S) 

ANCILLARY    SUEROUTIN>r 
SUBTRACT    A    SCALAR     FROM 
X(N2,I )=X(N1,I)-S 

REAL*£    X(20,101),S 

COMMON    XtN 

00    1     I=L»N 

X(N2t  I)=X(N1,I )-S 

RETURN 

END 


EACH    TERM    OF    A    VECTOR 


C 
C 
C 


SUBROUTINE    MULS  (Nl , N2 , S ) 
ANCILLARY    SUBROUTINE 
MULTIPLY    EACH    TERM    OF 
X(N2,I )=X(N1  ,1 )*S 

REAL*8    X(20t101),S 

COMMON    X»N 

00    1     1=1, N 

X(N2,I )=X(N1,I >*S 

RETURN 

END 

SLBROLTINE    DI VS (N1»N2 fS ) 
ANCILLARY    SUBROUTINE 
DIVIDE    EACH    TERM    OF    A 
X(N2,I  )  =  X(N1  ,1  )/S 

REAL*8    X(20fl01)vS 

CONMQIS    X,N 

DG    1    1=1, N 

X(N2, I)  =  X(N1,  I  ) /S 

RETURN 

END 

SUBRCUTINE    DUPV(N1,N2) 
ANCILLARY    SUBROUTINE 
DUPLICATE    A    VECTOR 
X(N2,I )  =  X(N1  ,1 ) 

REAL*8    X(20,101  ),S 

COMMON    X,N 

DO    1     1=1, N 

X(N2, I )=X(N1,I ) 

RETURN 

END 


A    VECTOR    8Y    A    SCALAR 


VECTOR    BY    A    SCALAR 


SLBRO 
AN 
IN 
REAL* 
CONMO 
EN=N- 
SN  =  1. 
R  =  E.M/ 
NINO  = 
NTNO  = 
FIVO  = 
THT3  = 
X(N2, 
X  (N2, 
NM2  =  N 
DO  1 
KP1=K 
KP2=K 
KP3  =  K 
ADD  = 
X(N2, 
X(N2, 
1  +  R*X ( 
RETUR 
ENC 


UTIN 

CILL 

TEGR 

8    X( 

N    X, 

1 

0  +  0/ 

2.4D 

R*9. 

R*l. 

R*5. 

R*l. 

1)=C 

2)=N 

-3 

K  =  l, 

+  1 

+  2 

+  3 

THTO 

KP2) 

N)=X 

N1,N 

N 


E     INTV(N1,N2) 

ARY    SUBROUTINF 

ATE    A    VECTOR    USING    MILNE'S    METHOD 

20,101 ), EM, R ,ADD,NIN3,NTN0,FI70,THT0 

N 

E'J 

+  1 

D  +  0 

9D  +  1 

D+0 

3D  +  1 

.0  +  0 

IN0*X(N1,1  }  +  NTN0*X(Nl,  2)-FIVC*X(Nl,3)+R*X(Nl,' 


NM3 


*(X(N1,KP1)+X(N1,KP2))-R*(X(N1,K)+X(N1,KP3)) 

=X<N2,KP1 )+ACD 

(N2tN-l)+NIN0*X(NlfN)  +  NTN0*X(M  ,N-1)  -FIVC*X(N1 

-3) 
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Appendix  C 
Examples 

The  listing  on  the  next  page  is  of  a  MAIN  program  which 
calls  RODISK  twice  to  solve  two  different  problems.  On  an  IBM 
360/67  compile  time  (for  the  MAIN,  RODISK,  and  the  ancillary 
subroutines)  was  12  s,  link  time  was  2  s,  and  execution  time 
for  both  broblems  was  1.5  s,  In  both  problems  we  used  N  =  hi 
and  iterated  10  times. 

The  first  problem  is  of  type  2  (solid  disk)  with  b  =  10, 
t  =  l/(1.6+.008r2),  v  =  0.3,  W   =  120,  a  (b)  =  14000,  and 
EaT  =  25105  +  1300  log  (t)  -  r2(233+l6t)  Units  are  inches  and 
pounds.  The  problem  was  made  up  from  the  exact  solution 

a  =  9000  +  50r2;    aQ  =  a     +  r2(120+l6t) 

The  results,  shown  on  page  18,  show  evaluations  for  the 

stress  components  which  are  correct  within  0.03  psi  even  though 

convergence  was  complete  to  only  about  4  digits  as  indicated  by 

the  sequence  of  values  above  the  final  tabulation;  these  are 

values  of  a   (0)/a  (b)  and  of  B. 
r    r 

The  second  problem  is  of  type  1  with  a  =  2.57,  b  =  5.15, 

a   (a)  =  18205,  a   (b)  =  22000,  t  =  0.1^93  r    ,  T  =  60  -  1.6r2, 

v  =  0.3,  and  Ea  =  19*1.3.  The  units  are  inches  and  pounds. 

Since  the  thickness  variation  is  a  power  law,  the  theoretical 

solution  given  in  the  body  hereof  may  be  used  to  obtain  the 

theoretically  exact  solution 

a    =  -113. 95r2  +  15832^  -  11708rQ 
r 

aa   =  122. 80r2  +  13801^  +  l4310rq 
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2  0 


T 
E 
K 
A 
3 
S 
S 
8 
0 

Y 
X 
X 
R 
X 
T 
X 
X 
X 
30  C 

c 

s 

E 


J. 

TH 
MPLI 
EAL* 
NTEG 
GMMO 
CMMO 

TH 
P(l) 
P(2) 
P(3) 
P(4) 
RAD  = 
RB-1 
=  41 
NM  =  N 
0  20 
IM=I 
=  EIM 
HICK 
(2), 
(2,1 
=  2.5 
(3,1 
CNTI 
=  X(2 
EEBE 
ALL 
ALL 

TH 
EEBE 
NM  =  N 
P(l) 
RAD  = 
RAO= 
RA=1 
RB  =  2 
ETA  = 
0  30 
IM=I 
=  5IN 
(5,1 
(6,  I 
=  ARA 
(4,1 
HICK 
(1,1 
(2,1 
(3,1 
ONTI 
ALL 
TCP 
NO 


E.  BR 

IS  IS 
CIT  RE 
3  X(2J 
ER  KP( 
N  X,N 
N  /ONE 
5  FIRS 
=  2 
=4 
=  10 
=  1 

1.0  +  1 
.40+4 

-1 

1=1, N 
-1 

*8RAD/ 
=1.0+0 
I  )=THI 
)=120. 
1050+4 
)  =  W/SR 
NUE 
C,N) 
E  =  S 

CIVS(2 
RODISK 
E  SECQ 
E  =  1.49 
-1 
=  1 

2.570+ 
5.150+ 
.82050 
.20  +  4 
5.D247 

1=1  ,N 
-1 

/EMM 
)=Y 
)=Y 

D+(BRA 
)=R/BR 
=1.493 
)=THIC 
)=BETA 
)=(6.0 
NUE 
ROOISK 


OCK  1    ^AY    1978 

A    MAIN    PROGRAM    TO 
AL*8(A-H,0-Z) 
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where  p  =  .29171  and  q  =  -1.87171.  The  computer  results  agree 
with  five  digit  accuracy  even  though  the  sequence  of  iterants 
(A,B)  shows  only  four  digit  convergence. 

The  same  two  problems  were  also  worked  with  N  =  101  and  14 
iterations.  The  execution  time  went  from  1.5  s  to  3.0  s.  The 
maximum  change  in  any  stress  value  was  0.3  psi.  Prom  these  and 
other  problems  it  may  be  concluded  that  there  are  no  difficulties 
of  accuracy,  computer  storage,  or  execution  time. 
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